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Given the sensitivity of current ground-based Gravitational Wave (GW) detectors, any continuous- 
wave signal we can realistically expect will be at a level or below the background noise. Hence, any 
data analysis of detector data will need to rely on statistical techniques to separate the signal from 
the noise. While with the current sensitivity of our detectors we do not expect to detect any true 
GW signals in our data, we can still set upper limits (UL) on their amplitude. These upper limits, 
in fact, tell us how weak a signal strength we would detect. In setting upper limit using two popular 
method, Bayesian and Frequentist, there is always the question of a realistic results. In this paper, 
we try to give an estimate of how realistically we can set the upper limit using the above mentioned 
methods. And if any, which one is preferred for our future data analysis work. 
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I. INTRODUCTION 

Gravitational Waves (GWs), ripples in space-time 
which travel at the speed of light, are a fundamental con- 
sequence of Einstein's General Theory of Relativity. Due 
to the great distance to any likely detectable sources of 
GWs, the signal amplitude reaching us will be very small. 
Because of limitations in technology, there have been no 
direct detections of GWs so far. However, with the future 
generation of detectors, like Advanced LIGO, we should 
be able to detect a variety of sources. 

Because of their nature, continuous GWs (emitting 
from axisymetric rotating neutron stars) reaching Earth 
are expected to be extremely weak. Therefore even with 
the quite significant sensitivity of our current detectors, 
it will be difficult to detect them. One possible way to 
increase the overall signal compared to the background 
noise (signal-to-noise ratio (SNR)) is to coherently inte- 
grate the data for several days up to few years. 

The basic problem in GW detection is to identify a 
gravitational waveform in a noisy background. Because 
all data streams contain random noise, the data are just 
a series of random values and therefore the detection of 
a signal is always a decision based on probabilities. The 
aim of detection theory is therefore to assess this proba- 
bility. 

The basic idea behind the current methods of signal 
detection is that the presence of a signal will change the 
statistical characterization of the data x(t), in particular 
its probability distribution function (pdf) P(x). Recall 
that the pdf is defined so that the probability of a random 
variable Xi lies in an interval between x(t) and x(t) + 
dx is P(x)dx. Let us denote by P{x\Q) the probability 
of a random process x{t) (representing our data) in the 



absence of any signal, and by P{x\h) the probability of 
that same process when a signal h(t) is present. Given a 
particular measurement x(t) obtained with our detector, 
is its probability distribution given by P(x\0) or P(x\h)l 
In order to make that decision, we need to make a rule 
called a statistical test. 

There are several approaches to find an appropri- 
ate test, notably the Bayesian, Minimax and Neyman- 
Pearson approach (for an overview, we refer the reader 
to Jaranowski and Krolak, |l| and the references listed 
therein). In the end, however, these three approaches 
lead to the same test, namely the likelihood ratio test 

Among the three main approaches, the Neyman- 
Pearson approach is often used in the detection of gravi- 
tational waves 3 • This approach is based on maximizing 
the detection probability (equivalently minimizing the 
false dismissal rate) for fixed false alarm rate, where the 
detection probability is the probability that the random 
value of a process which contains the signal will pass our 
test, while the false alarm probability is the probability 
that data containing no signal will pass the test nonethe- 
less. Mathematically, we can express the Detection and 
False Alarm probabilities as Q 



(1) 
(2) 



P D {R) = / P(x\h)dx, 
Jr 

P F {R) = [ P(x\0)dx, 
Jr 



respectively, where R is the detection region (to be de- 
termined) . 

The Likelihood ratio A is the ratio of the pdf when the 
signal is present to the pdf when it is absent: 



A = 



P«t)\h(t)) 
P(x(t)\0) 



(3) 



* Email: iraj.gholami@theorie.physik.uni-goettingen.de 



Taking the data to be x(t) = h(t) + n(t), with h(t) 
the signal and n(t) the noise and with the assumption 
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that the noise is a zero-mean, stationary and Gaussian 
random process, we can write the likelihood ratio as 



A = 



P(x\h) 
P{x\0) 

cxp(— i(x — h\x — h)) 
exp(—^(x\x)) 

exp[(x\h) - hh\h)]. 



This leads to the log of likelihood function as 
logA = (x|/i)-i(/i|/i). 



(4) 



(5) 



We can also rewrite the simple expression of Eq. ([5]) 
for the likelihood function in terms of the new variables, 
A a and h„ 



as 



1 



log A = (x\A a h a )--(A a h a \A"h b ). 

where the constant (in time) amplitudes A a 
A a (h ,ip,i,$> ) are H 



.4 1 
A 2 
A 3 
A 4 



A + cos $o cos 20 — A x sin $o sin 20, 
A + cos $o sin 20 + A x sin $o cos 20, 
— A + sin$o cos 20 — A x cos$o sin 2-0, 



(6) 



(7) 
(8) 
(9) 



= —A + sin$ sin 20 + A x cos $ cos 20. (10) 



and 



hi(t) 
h 2 {t) 
h 3 (t) 
hi{t) 



a(t) cos0(t), 
b(t) cos (j>(t), 
a(t) sin 4>{t), 
b[t) sin </>(*), 



(11) 
(12) 

(13) 
(14) 



where a(t) and b(t) are functions of right ascension a 
and declination 5; they are independent of 0, and is 
the phase of the wave signal seen at the Solar System 
Barycenter (SSB) |5|. Likewise 



A+ = -h (l + cos 2 l), 
A x = ho cos i, 



(15) 
(16) 



where ho is the wave amplitude, i the inclination angle, 
the polarization angle and $o the initial phase. 

Since the A a s depend neither on the detector proper- 
ties nor on the frequency or the time, we can take them 
out of the inner product and write the log of likelihood 
ratio as 



log A = A a (x\h a )-^A a A b (h a \h b ). 



Defining the new variables 



(x\h a ), 



(17) 



(18) 



and 



we have 



M ab = (h a \h b ), 



\ogk = A a H a --A a A b M ab . 



(19) 



(20) 



The maximum detection probability follows from the 
maximization of the likelihood function: by maximizing 
the likelihood function with respect to the A a (which, 
again, are independent of the detector), we have 



d log A 

dA a 



0. 



This leads us to 



and therefore 



Hn. A,. T „M n 



< LE = (M~ 1 ) ab H a 



(21) 



(22) 



(23) 



The label AILE denotes the Maximum Likelihood Esti- 
mator; it corresponds to the values for the A a s we cal- 
culate from our data by maximizing the likelihood ratio 
(so that, in practice, we are calculating A a — E[A a MLE \). 
By definition, the ^-Statistic is the maximum of the log- 
arithm of likelihood function. Substituting Eq. (|23p into 
Eq. (I20l). we have 



T = log A 



~ H a (M- V ) ab H b 



(24) 



This is the ^-Statistic which a generalized version of that 
for the multi-IFO (Interferometer Observatory) can be 
found in Q. 

Again by writing the data as x(t) = n(t) + h(t), and 
using Eq. (2.5) of Cutler-Schutz (CS) [4] indicating that 
((x\n)(y\n)) — (x\y) with this fact that ((h\n)) — 0, we 
would have the following 



(2T)=4+(h\h), 



(25) 



which follows a x 2 distribution with 4 degrees of free- 
dom and non-centrality parameter p 2 = (h\h), that is 
the square of optimal signal to noise ratio (SNR 2 ). The 
degrees of freedom come from the 4 unknown parameters 
of pulsar namely the amplitude (ho), inclination an- 
gle (i), polarization angle (ip) and the initial phase 
($o)- As was described in [J], even for a multi-IFO the 
number of freedom will remain unchanged, since we al- 
ways have the same 4 unknown parameters in entire the 
search. 

The ^-statistic shown above, which was originally de- 
rived by Jaranowski, Krolak and Schutz (JKS)[5[, is the 
optimal statistic for detection of nearly periodic gravita- 
tional waves from GW pulsars. We can use this statistical 
tool to search for any kind of pulsars; unknown sources 



3 



or targeted search. In targeted search we know every- 
thing about the source by using some other astronomical 
techniques, such as radio astronomy, gamma and X-ray 
astronomy etc. The main information required for our 
search are; frequency and its derivatives and position of 
the source. Given all these information to the software 
developed by the LSC (LIGO Scientific Collaboration) 
|7J (the implementation to our work can be found in 
), and using a single workstation, we will be able to 
search a known pulsar in a few minutes. 

The work in this paper was done using simulated data 
for 100 arbitrary pulsars. To set the positions and the 
frequencies of these simulated pulsars, we have used the 
data of the known pulsars, as given in the Australian 
Telescope National Facilities (ATNF) catalogue @. To 
make the simulated data more realistic, we generated the 
data at the level of LIGO detectors sensitivity. Since 
our detectors (initial LIGO and even Enhanced LIGO) 
are not sensitive enough to detect any signal until now, 
we assume that the data are just simply noise without 
any signal in them. Due to that, our simulated data are 
just pure noise and therefore instead of looking for any 
detection, we set the upper limit on the strength of the 
gravitational wave signal. For the historical reason we 
take the value of 95% for the upper limit. All the search 
done in this paper are in frequency domain. As the main 
goal of this paper is to compare two different approaches 
in setting upper limits, we use Bayesian and Frequentist 
algorithm to perform it on the same data for each pulsar. 
These two algorithms will be explained in more details 
below. 



II. FREQUENTIST UPPER LIMITS 

The frequentist probability of an event represents the 
expected frequency of occurrence of that event. The re- 
sult for our upper limits depends crucially on the exper- 
imental data under examination. The confidence value 
associated with these upper limits indicates the expected 
occurrence of detection statistics values more significant 
than the one that we have measured in the presence of 
signals whose amplitude is equal to the upper limit value. 

To set the frequentist upper limit on the amplitude of 
gravitational waves, we use the .F-statistic as an optimal 
detection statistic. To start with, we need to assign a 
confidence level C - roughly speaking, our criterion will 
then be that, for our "repeated measurements", in C- 
percent of the time the value of 2T is above a specified 
threshold. 

Let us explain how this works in detail for the example 
of setting a 95% upper limit on ho- For this, we need to 
find at which ho it is true that 95% of the values of the 
.F-statistic are above the initial value of 2.F derived from 
the data. To do so, we proceed step by step as follows: 

1. Compute the .F-statistic of a perfectly matched sig- 
nal using the exact values for the signal parame- 



ters (such as frequency, longitude, latitude and fre- 
quency derivatives). Let us call the resulting value 
of the .F-statistic T* . 

2. Estimate the signal amplitude, ho, using our param- 

eter estimation routine. 

3. Take this ho as the initial value of the search. 

4. Since we assume that there is no signal in the data - 

that it is pure noise -, we can randomly assign ar- 
bitrary values to the other signal parameters (such 
as 4>o, ip and cos/,). 

5. To determine the probability distribution of T- 
statistic, we take a random frequency value with a 
band of 0.1 Hz around the actual pulsar frequency 
(as was proposed in the LIGO SI paper [l(|) and 
inject the artificial signal. With this choice, we are 
sure to be on the safe side; we use a large amount 
of data (in order of several months up to a year), so 
that the 0.1 Hz band will not lead to any spurious 
correlations between the search parameters. 

6. After injection, compute the .F-statistic once again. 

Let us designate the resulting value of 2.F as J 7 '; 
store this value for later use. 

7. Repeat the injections, computing of .F-statistic for 

150 times. Save all resulting values of J 7 '. (The 
number of iterations used here is a heuristic value.) 

8. As we are looking for a 95% upper limit, proceed as 

follows: if the confidence level (the percentage of 
instances in which T' is greater than J 7 *) was less 
than 90% or above 98% (say x), multiply the ho by 
the ratio of — and take this value as the initial ho 
for the next step. 

9. Repeat steps "6" and "7" until the confidence level 

is in one of the following ranges: a) 90% — 95%, or 
b) 95% - 98%. 

10. For case a), multiply ho by 1.05; for case b), multi- 

ply by 0.90. (The factors 1.05 and 0.90 are, again, 
heuristic.) 

11. Repeat the calculations of step "7" and following, 
but this time with 1000 injections in each run (in- 
stead of 150) to improve the statistics. 

12. Repeat step "11" for 6 times; in each run follow 
the instructions in step "10" . (The number of rep- 
etitions is heuristic; it is chosen in a way that the 
range of computed confidence levels will always in- 
clude values higher and lower then 95%; therefore 
we can make an " interpolation" fit instead of hav- 
ing to extrapolate.) 

A flow chart version of this procedure can be found in 

Fig. m 
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Compute the 2F value of a 
perfectly mtached signal (F*) 



Estimate the value of hO and 
take that as an initial value 
to search 



Inject the signal randomly 
distributed in frequency band - 
of 0.1 Hz 



Compute the 2F value (F') 



Save the value in a file 



Multiply hO by 95/X 



No 

h 



r X = N(F'>F')/ N 
90 < X < 98 



Multiply hO by 1.05 



Inject the signal randomly 

j »- distributed in frequency band 

T of 0.1 Hz 



Multiply hO by 0.95 



Repeat 
1000 times 



Inject the signal randomly 
distributed in frequency band 
of 0.1 Hz 



Compute the 2F value (F) 
and store it 



Compute the 2F value (F) 
and store it 



Repeat 
6 times 



Compute X 
and store it 



Compute X 
and store it 




Interpolate the 6 final X 
to get theh0(95%) 



FIG. 1: A Flow-Chart of how we implemented the frequentist 
upper limit. 



III. B AYE SI AN UPPER LIMIT 

The Bayesian probability is a measure of degree of belief 
in the occurrence of a statistical process. In contrast with 
the Frequentist probability, in the Bayesian approach, we 
do not need for an event of that particular type to have 
actually happened; all we need is to find a measure for the 
degree to which a person believes that a given proposition 
is true. 



Our goal is to set an upper limit on the strength of the 
gravitational wave amplitude, /iq, using a given amount 
of available data - which is the posterior probability of 
ho that we look for. Therefore, our data plays the role of 
B in Eq. (f^C?)) ; which we denote it by s. The term A is 
the quantity about which we intend to draw conclusions 
using our data; in our case, this is the upper limit /io, so 
we will substitute ho for A in what follows. With these 
substitutions, Eq. J2H]) now reads 



P(ho\s) 



P(s\h Q ) X P(hp) 
P(s) 



(27) 



where P(ho\s) is the conditional probability of ho (pos- 
terior probability) given the data s, P(s\ho) is the like- 
lihood function (to be defined below), and P{ho) repre- 
sents our prior knowledge about the distribution of ho- 

Since the term P(s) is independent of our signal, we 
can consider it as the constant normalization factor; it 
will cancel out automatically when we compute the con- 
fidence level. Therefore, we can rewrite our Bayes' theo- 
rem for the general case of all signal parameters as 

P(h , 4>, h $o|s) oc P(s\h , ip, i, $o) x P(ho, t, $o), 

(28) 

where, again, P(ho,il>, t, $o| s ) is our posterior probabil- 
ity (to be calculated), P(s\ho, tp> L > ^o) is the likelihood 
function and P{ho, ip, i, $o) is the prior probability of 
h ,ip, t,®o- 

There are two common choices for estimating prior 
probability, known as a flat prior and Jeffrey's prior. 
In the flat prior, the prior probability is chosen to be 
constant (P(ho) = constant), while in Jeffrey's prior, it 
is taken to vary inversely proportional to the value of ho 
(P(h ) = I /ho). For more details we refer the reader to 
[T3| and a comparison for this case can be found in [l3j]. 

The Jeffrey's prior gives a higher value in upper limit 
than the flat prior, while a flat prior gives a more realistic 
value for our case [HI . Therefore in the following we will 
focus on a flat prior, as the case followed in 

To obtain the posterior probability, we need to cal- 
culate the likelihood function. By Eq. (IM|) . it can be 
expressed as 



A. Theoretical approaches 



P(s\ho,*P,L,<Z> ) oce-i M ^ Aa - Aa °X Ab - Ab °) =G, (29) 



The key ingredient of the Bayesian approach is the 
Bayes' theorem (a simple proof of that can be found in 



P{A\B) 



P{B\A)P(A) 
PjB) ' 



(26) 



The term at the left hand side is called the posterior prob- 
ability, while P{A) is the prior probability which reflects 
our initial knowledge about the quantity A. The term 
P(B\A) is called the likelihood function; the log of this 
is, in fact, the J r -statistic to be computed from our data. 



where A = (A 1 , A 2 , A 3 , A 4 )(h , ip, h $o) are the four 
amplitude parameters defined in Eqs. ([T UTU)) and G = 
G(ho,if>,i,$o). The A = (A lf > 7 A 2 ", A 3 °, A 4 °) are also 
the best fit for the A a s resulting from our calculation of 
the ^-statistic. 

For proper normalization, we first compute the integral 

/>oo />1 / >7r /4 p2tt 

1=1 P{ho)dho da di> / d$ G, (30) 

JO J-l J-tt/4 Jo 

where /i = cosi, and we will set l P(ho) = constant'. To 
find the upper limit we use h r nax as the upper bound in 
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the integration over ho, 

Iul = / P(h )dh / dfi / # / d$ G. 

Jo J-i J-tt/4 Jo 

(31) 

We select /i™ oa: in such a way that the ratio Iul /I gives 
us the desired confidence level. In our case, we are look- 
ing for the 95% upper limit, therefore 



Iul 
I 



0.95. 



(32) 



B. Practical implementation 



To implement the above formalism, let us first con- 
struct the function G(h , ip, & o)- To do so, we need 
to expand the matrix of Eq. (|19l) . The elements of this 
matrix depend on the three amplitude modulation coeffi- 
cients (A, B and G) defined in [J]. Based on the notation 
used here, these elements take the form of 



A/2 C/2 
C/2 B/2 

A/2 C/2 

C/2 B/2 



(33) 



which a detailed procedure of their derivation can be 
found 81. Then we can construct the four elements 



Gi = 


4<* 


-A 1 ") 2 -\ 




- A l0 )(^ 2 


-A 2 °),(34) 


G 2 - 




- A 2 ") 2 - 




-A 2 °)(A X 


-A 1 ") (35) 


G 3 = 




- A 3a ) 2 -\ 




-A 3o )(yl 4 


- A 4 °),(36) 


G 4 = 


> 


-A 4 ") 2 - 




- A 4o )(A 3 


- A 3o ){37) 



to make the final form of G(ho, ip, t, $o) m Eqs. (|3"0"|) and 
J2B as 



Bayesian UL 
Frequentist UL 



5 10 15 20 25 30 35 40 45 50 55 60 65 70 75 80 85 90 95 100 

Pulsar Numbers 



FIG. 2: The upper limits value on the ho for 100 arbitrary 
pulsars using simulated data. In this plot, both Bayesian and 
Frequentist ULs are shown. 



satisfied. This can be done using the Numerical Integra- 
tion routines in mathematical software like Mathematica 
and Maple. Another way would be to calculate the pos- 
terior probability of ho by marginalizing over the other 
three parameters. This can be expressed in mathematical 
form as 



p(h \s) oc 



G(h Q , ip, i, $ ) dip dfi d$ - (39) 



Once the posterior probability for ho is known, one can 
then integrate it over a sufficient range of ho to find out 
the area covered; the result can be used for proper nor- 
malization (namely unit total area). Next, we can find 
out at which ho the fraction of area would satisfy our 
required confidence level. 

Both the above methods have given equivalent results 
as discussed in details in [8j. However, for the work ex- 
pressed in this paper, we followed the second algorithm. 



IV. RESULTS AND DISCUSSION 



G = cxp[--(G 1 + G 2 +G 3 + G 4 )]. (38) 

This is the core equation for our upper-limit analysis in 
Bayesian approach. To construct this, we need all the 
above mentioned parameters to be resulted from our soft- 
ware. The software we have used for this purpose was 
developed partly by the author of this paper and is now 
part of the LAL (LIGO Algorithm Library) @. With this 
software we calculate the four amplitudes A a as well as 
the matrix elements M a b (namely the amplitude modula- 
tion coefficients A, B and C). Once we have constructed 
the likelihood function GQiq, ip, l, $o), we can calculate 
the UL value in Eq. (|3"2")l in two ways. One is to follow 
the exact procedure spelled out above; first calculating 
the normalization in Eq. (|30[) and then trying to find 
a value of h™ ax for which the ratio of Eq. (|3"2"j) will be 



Once again, we have selected 100 arbitrary pulsars fre- 
quencies and positions (based on the real pulsars infor- 
mation taken from ATNF [j|). We have generated the 
simulated data at the level of current LIGO detectors 
sensitivity (using the LAL software @, 0] developed by 
LSC) and computed the UL for these pulsars. To com- 
pare with the real data, the simulated data contains just 
noises where the upper limit set on them are shown in 
Fig.H 

The blue rectangular in this plot represent the value 
of upper limits in Bayesian approach and the red circle 
points to the Frequentist ones. The horizontal axis indi- 
cates the pulsars number, therefore on each vertical line 
corresponds to each pulsar we should have one blue rect- 
angular and one red circle. However, as seen, in some 
cases there is just one blue rectangular and missing red 
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FIG. 3: One of the bad upper limit value on the ho in Fre- 
quentist approach. The value of IT for this case was 0.08. 



circle (Frequentist UL). These count for 16 pulsars, in 
which 7 of them were caused due to some unknown rea- 
son. The reason for the 9 others is that; in these cases, 
due to very small amount of 2T, the Frequentist upper 
limit procedure could not be converged to any particular 
value. It means, for some cases, by changing (increasing 
or decreasing) the value of ho for more than two order 
of magnitude, the upper limit value always stands above 
96% or 97%. 

An example of that is shown in Fig. [3j that is the case 
where 2T = 0.08. This figure shows the dependency of 
Frequentist upper limit to ho for an individual pulsar. 
It's clear that even by increasing the ho for about 2.5 
order of magnitude, the value of upper limit lies mostly 
about 100%. While, in general the upper limit is very 
sensitive to small changes in ho- This was done 76 times 
and in each time the value of ho was increased according 
to the procedure expressed in Sec. HU Note that the 
starting and ending value of ho are significantly smaller 
than the ho require for 95% upper limit shown in Fig. [5] 
Means that, in normal condition where the required ho 
to get 95% upper limit is in order of 10 -26 , by setting the 
ho in the range of 10 -28 — 10 -27 , we should get a very 
small upper limit compare to 95%. In fact, as disscused 
below, the low value of 2T for this pulsar is the reason 
of such a behavior in the Frequentist framework. 

The same behavior was shown in [H, [l4[ by using the 
real data. The reason is clear; we have pointed out that 
the Frequentist approach is based on the number of occu- 
rance of an event. For this we should set a threshold and 
count how many times the value of that particular pa- 
rameter is passing this threshold. Naturally there can be 
some False Alarms (FA) , which a noise shows itself strong 
enough to pass this threshold. Since the ^-statistic fol- 
lows a x 2 distribution with four degrees of freedom, the 
FA follows as (equation 3.44 of @), 

a= (1 + 2T) c- 2jr . (40) 

The values of 2T in which the Frequentist upper 
limit could not be converged are: 0.08,0.34, 0.46,0.47, 
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FIG. 4: The ratio of Bayesian over Frequentist upper limits 
value on the ho for 100 arbitrary pulsars using simulated data. 



0.49, 0.51, 0.59, 0.63 and 0.84. So, by putting these num- 
bers in the above equation we get a very large values 
of FA. For example, in the case of 2T = 0.84 we have 
a = 0.80, 2T = 0.46 gives a = 0.92 and for 2T = 0.08 
we get a = 0.997 ~ 100. This explains the whole story. 

This high false alarm probability means that noise 
alone has a high chance of producing an J-"-statistic value 
greater than the J-* produced by our data set and our 
template. In other words, our realization of the noise is 
one that; it is particularly unlikely to look like it con- 
tains our signal, and this statistical fluctuation yields a 
low upper limit value or even does not converge - due 
to the nature of the procedure that we use to determine 
the upper limit. The Bayesian approach is less sensitive 
to such fluctuations. Note that the resulting frequentist 
upper limit is not an artifact of our technique, but it is 
still a perfectly correct and consistent upper limit in the 
Frequentist framework. This clearly shows the nature of 
our data. 

To investigate further and check the behavior of the 
Frequentist UL, let us compare its value with that of 
Bayesian. To do so, we plot the ratio of the Bayesian UL 
over the corresponding value of Frequentist for the same 
pulsar versus the value of 2T (Fig. H]). The output says; 
as we go further to the value of 2T less than 4, the ratio 
increases and when we reach to the 2T = 1, this ratio is 
quite significant; about a factor of 2.3. For the case of 
2 J- < 1, we have already seen that the upper limit for the 
Frequentist approach did not converge. Therefore they 
are not shown in this plot. If they would converge, they 
should be quite lower than that of Bayesian and therefore 
the ratio should go much higher. The same behavior was 
shown in Q with real data. 

Fig. [4] also shows that, at roughly 2T — 4 the ra- 
tio is close to unity and roughly remains the same when 
2T > 4. This tells that the problem of low value in 
upper limit in Frequentist approach appears when we 
have 2T < 4 while for larger value of .F-statistic there 
is always agreement between Frequentist and Bayesian 
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frameworks. 

Apart from the difference in the nature of Frequentist 
and Bayesian frameworks, there is another difference in 
performing a Frequentist and the Bayesian upper limit 
search. Since in Frequentist algorithm we need to in- 
ject some artificial signals into the data and then search 
the newly generated data to compute the ^-statistic in 
each iteration, this requires a high amount of computa- 
tional resources. To increase the sensitivity we need to 
use more data that requires more computational power 
as well. Because, the required time to search the data 
to compute the ^-statistic is linearly proportional to the 
amount of data. In order to have a better statistic in Fre- 
quentist algorithm, we therefore need a larger iteration. 
This would additionally brings another linear increment 
in the cost for the computation. While in the Bayesian 
approach, to compute the ^-statistic and the other com- 
ponents, we search the data just once. Then compute 
the P(ho) by marginalizing the probability over the if), i 
and $o- These all will be done once and are computa- 
tionally very cheap. As an estimation, the entire process 
for one pulsar using Bayesian algorithm takes about half 
an hour up to one hour in a single workstation. In a 
good approximation this is independent of the amount 
of data. Because, searching in the large amount of data 
(say about one year) to compute the ^-statistic and other 
components takes just about few minutes. In contrary, 
the required time for a Frequentist algorithm to search 
for single pulsar in an amount of data in order of one 
year takes about 3 weeks on a single workstation. 

As a summary; although search in the Frequentist up- 
per limit shows the exact nature of our data, however 
there are some disadvantages with the same search by 
using the Bayesian algorithm. The important one is that 
in the case where our data shows a small value of T- 
statistic in a particular frequency bin and position of the 
pulsar, we cannot trust the upper limit value produced 
by Frequentist approach. Likewise, performing a search 
in the Frequentist framework is much expensive than the 
same search in Bayesian approach. 



V. DISTRIBUTION OF 95% BAYESIAN UPPER 
LIMITS ON h USING SIMULATED DATA IN 
FREQUENCY DOMAIN 



95% upper limit on ho of each data set. The sky loca- 
tions in the search are chosen randomly such that their 
distribution over the solid angle is uniform; detectors po- 
sition are picked randomly from a list consisting of the 
locations of HI, LI, VIRGO and GEO600 detectors. The 
resulting mean upper limit is 



95% \ 



(10.67 ±0.04) 



s h (f) 



T 



(41) 



In order to compare our result with the simulation in 
Time Domain (TD) done by Dupuis and Woan |15| . we 
repeated this experiment with only HI, LI and GEO600 
detectors, as in their analysis. The results are in a very 
good agreement with a ratio in ULs 



K % ) 



FD 



= 0.98. 



(42) 



TD 



Fig. [5] shows the distribution of for 5500 differ- 
ent runs in frequency domain, which is also in a good 
agreement with that of presented in [15j. 

These results show that although we use different do- 
main (TD or FD) to search for gravitational waves, if 
we stay in Bayesian framework, both give the same 
results theoretically (a more details can be found in 
8]). However, using different frameworks (Frequentist 
or Bayesian) will may lead to a different outcome. 



n 



20 25 



As an application of Bayesian algorithm, we now 
present the distribution of ULs (95% upper limit on ho) 
computed in the Frequency Domain (FD). We start with 
the idealized case of a large number (5500) of simulated 
data set with pure noises (no signal), and compute the 



FIG. 5: Distribution of 95% Bayesian upper limit on ho us- 
ing 5500 individual simulated data runs in frequency domain. 
(hf % ) = (10.59 ±0.04). 



[1] P. Jaranowski, A. Krolak, Living Rev. Relativity 
8,(2005), 3. 

[2] M.H.A. Davis, "A review of Statistical Theory of Signal 
Detection", in B.F. Schutz, ed., Gravitational Wave Data 



Analysis, Proceeding of the NATO Advanced Research 
Workshop, held at Dyffryn House, St. Nicholas, Cardiff, 
Wales, 6-9 July 1987, vol. 253 of NATO ASI Series C, 
73-94. 



8 



[3] Bernard, F. Schutz, Introduction to the Analysis of Low- 
Frequency Gravitational Wave Data, gr-qc/9710080 

[4] Curt Cutler and Bernard F. Schutz, Phys. Rev. D72, 
063006 (2005). 

[5] P. Jaranowski, A. Krolak, and B.F. Schutz (JKS), Phys. 

Rev. D58, 063001 (1998). 
[6] The LAL software can be found on the following websites: 

http : //www. lsc-group .phys .uwm. edu/daswg/ 

projects/lalapps .html 
[7] The doxygen version of the LAL softwares can be found 

in this site: 

http : //www. lsc-group .phys .uwm. edu/lal/ 
slug/nightly/ doxygen/html/ dirs . html 

[8] Iraj Gholami, PhD thesis, Potsdam University, Germany 
(2008) , |http://opus.kobv.de/ubp/volltexte/2008/1880/ 



[9] The Australian Telescope National Facility catalogue 

http://www.atnf.csiro.au/research/pulsar/psrcat/ 
[10] B. Abbott et al. (The LIGO Scientific Collaboration), 

Phys. Rev. D69, 082004 (2004). 
[11] Glen Cowan, Statistical Data Analysis, Oxford Science 

Publications, 1998. 
[12] P.C. Gregory, Bayesian Logical Data Analysis for the 

Physical Sciences, Cambridge University Press, 2005. 
[13] Rejean Dupuis, PhD Thesis, University of Glasgow 

(2004). 

[14] Iraj Gholami, PhD thesis defense presentation, 2008. 
[15] Rejean Dupuis and Graham Woan, Bayesian estimation 

of pulsar parameters from gravitational wave data, Phys. 

Rev. D72, 102002 (2005). 



